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ABSTRACT 


The multigrid method (MGM) has been used to numerically solve the 
pair of nonlinear elliptic equations commonly used to generate two- 
dimensional boundary-fitted coordinate systems. Two different geometries 
are considered: one involving a coordinate system fitted about a circle and 
the other selected for an impinging jet flow problem. MGM uses a nest of 
grids from finest (upon which the solution is sought) to coarsest and is 
based on the idea of using relaxation sweeps to smooth the error (equivalent 
to eliminating high frequency Fourier components of the error) . Thus most 
of the computational work is done on coarser subgrids to eliminate longer 
wave length components of the error. Two different relaxation schemes are 
tried: one is successive point overrelaxation and the other is a four-color 
scheme vectorizeable to take advantage of a parallel processor computer for 
greater computational speed. Results using MGM are compared with those 
using SOR (doing successive overrelaxations with the corresponding relaxation 
scheme on the fine grid only) . It is found that MGM becomes significantly 
more effective than SOR as more accuracy is demanded and as more corrective 
grids, or more grid points, are used. For the accuracy required here, it is 
found that MGM is two to three times faster than SOR in computing time. With 
the four-color relaxation scheme as applied to the impinging jet problem the 
advantage of MGM over SOR is not as great. Perhaps this is due to the effect 
of a poor initial guess on MGM for this problem. 


The multigrid method (MGM) [1] can numerically solve linear or non- 
linear elliptic partial differential equations more rapidly than conven- 
tional means of solution such as successive overrelaxation (SOR) . MGM can 
be applied to the numerical solution of partial differential equations not 
amenable to numerical solution by fast direct matrix solvers such as 
diagonal decomposition. Thus it was deemed desirable to apply MGM to the 
numerical solution of the system of nonlinear elliptic equations commonly 
used to generate boundary-fitted coordinate systems, especially when the 
number of grid points is large. The standard elliptic equations for a typical 
mapping, shown schematically in Figure 1, are 

L (x,y) = ax - 2gx + yx + J 2 (Px + Qx ) = 0 (1) 
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and P and Q are functions of £ and n. Dirichlet conditions are specified on 
all boundaries of the computational space including the interior slit (which 
maps to the body in the physical space) . Each side of the slit has a set of 
Dirichlet data with a common value for each of the endpoints of the slit. 


The basic idea of MGM is to do most of the computational work on coarser 
corrective grids containing far fewer points than the finest grid upon which 
the solution is sought. The grids form a nest, each coarser grid having 
twice the mesh spacing in each coordinate direction of the previous finer 
grid. In Figure 2 which represents the Full Approximation Storage scheme of 
: u = (x,y), L = j such that Eqs . (1) and (2) become Lu = F = |q} ’ 


[ 1 ] 


1 = k = M (k representing the k th grid with M the finest), <t> = (x,y) on the 
boundaries of the computational space (Dirichlet values so that A is an 


[1] A. Brandt, Math, of Comp ., Vol. 31, No. 138, April 1977, pp. 333-390. 


identity operator) and superscripts refer to discretized quantities on the k 1 -* 1 
grid. (All operations involving 4) in the flow chart can be ignored, since the 
Dirichlet conditions are constant on all the grids.) The main idea behind MGM 
is that relaxation sweeps are a smoothing process which eliminate the highest 
frequency Fourier components of the error on any grid. First, starting with 
an initial guess for the solution, several sweeps are carried out on the 
finest grid to eliminate high frequency components of the error. The smoothed 
out error is represented by the residual f M = L M U^ 1 - L M u M and the correction 
U M -u M (where is the exact discrete solution on the finest M t ' 1 grid) . The 
residual, consisting mainly of longer wave-length Fourier components, is dealt 
with by solving its coarser-grid approximation 
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for U , which is represented by F for k = M-l in the lower right box of 
Figure 2. The symbol 1^ ^ means interpolation of a quantity from the k^ 
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grid to the (k-l)st grid. Eq. (4) is solved in the same way as the original 

equation on the finest grid. If solution of (4) is obtained after several 
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out function U -u is added to u . That is U -U + ( 
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is the expression in the lower left box for k = M. The new u is a better 
approximation to the solution and is the starting point for more relaxation 
sweeps for the original set of Eqs . (1) and (2) on the finest grid. If conver- 
gence is obtained, the process is complete; if not, the process returns to the 
coarser grid to sweep the residual equation again. If it doesn't converge 
after a few sweeps, then the next coarser grid is used to eliminate long wave 
length errors for the residual equation, etc. Each residual equation has a 
corresponding residual equation and correction on the next coarser grid. (The 
residuals were weighted locally as in [1].) 


Figure 3 shows computer drawn body-fitted coordinate systems generated 
to a specific accuracy using MGM and SOR (the two coordinate systems coincide). 
The relaxation scheme used was successive point overrelaxation. According to 
the notation used in Figure 1, m and n are 81 and 21, respectively; the slit 
end points are and ^53 ,r| 13^’ respectively; ( X £>Y b ) = ("8.4, -8.0) 
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and (x ,y ) = (7.6, 0.0); Ax = .2 and Ay = .4; A£ = An = 1; and the body is 
a circle of radius one centered at (x,y) = (0, -3.2). P and Q were set to 
zero in Eqs . (1) and (2). An experimentally determined, essentially optimum 
overrelaxation factor of 1.7 was used in the successive point overrelaxation 
sweeps in both MGM and the SOR method. All coarser corrective grids contain 
grid points on the slit. The initial guess for x(£,n), y(5,n) in the computa- 
tional space is obtained by extending the Dirichlet data at the outer boundaries 
throughout the space except at the slit, where the body Dirichlet data are used. 
The convergence criterion for the solution of Eqs. (1) and (2) was that both 

L 2 -error norms (one for each equation) be less than an input value j J E j | . 

" 2 

(This will be called satisfaction of ||e)| .) For Figure 3, j I E J I = .001. 
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To satisfy this criterion, MGM used 32.5 WU and 16.08 CP seconds compared to 
66.0 WU and 22,17 CP seconds for SOR. (A work unit (WU) is the equivalent of 
one SOR sweep on the finest grid, and CP seconds refer to central processor 
seconds used on the Texas Instruments Advanced Scientific Computer (TI-ASC) . ) 

For [ ] E | | = .01, MGM used 20 WU compared to 29 WU for SOR; CP time was the 

same for both methods (due mainly to the additional computational work in 
computing residuals in MGM) . The results show that the effectiveness of MGM 
increases (compared to SOR) as the error norm decreases. This is consistent 
with the fact that the remaining longer wave length errors are eliminated more 
slowly using SOR. The parameters 6 = .3, £ = .3 were used to control the flow 
of MGM. The parameter 6 determines the convergence test on each grid and the 
parameter C determines how fast the convergence must be (how fast the high 
frequency components are eliminated) on each grid. Whenever 


C < (||E||jM 1 + 1 /(||E! 
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(k-l)st grid with an error norm to be satisfied equal to 6 ( E ) 


(Super- 


scripts i,k refer to the ith relaxation sweep and the k c ^ grid, respectively.) 
These parameters are used as in [1], have a range (0 < 6 < 1 ; 0 < z, < 1) , and 
greatly influence the performance of MGM. The present choice is not necessarily 
optimum but was the best of a number of choices tried in the unit square. 

Figure 4 shows a computer drawn body-fitted coordinate system, similar 
to Figure 3, generated with MGM and satisfying ||e([ = .001. The grid 

parameters are (see Figure 1): m and n equal to 129 and 81, respectively; 
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slit end points of (^ g , n^) and (5 gl , n^ g ) , respectively; (x^,y b ) = 

(-7.68, -8.0) and (x^,y ) = (7.68, 0.0); Ax = .12 and Ay = .1; A£ = Ap = 1; 

and the circle of radius one was centered again at (0, -3.2). To satisfy 

] |e| | = .001 MGM used 21.863 ITU and 70.67 CP seconds compared to 102.0 WU 

l 2 

and 217.83 CP seconds used by SOR. This represents a significant saving of 

computer time by MGM. To satisfy | |eJ L = .01 MGM used 10.863 WU compared 

W 

to 17.0 WU used by SOR with CP time essentially the same. These results, 
along with those for Figure 3, show that MGM is more effective, compared to 
SOR, when more corrective grids are used and more accuracy is required. 

Figure 4 has five corrective grids and Figure 3 has three corrective grids 
(including the finest). The parameters 5 = .03 and 4 = .2 controlled MGM 
for Figure 4. Choosing smaller 6 and £ makes it more likely that all the 
coarser corrective grids will be used, which is desirable. 

Figure 5 shows a computer drawn body-fitted coordinate system generated 

using MGM and satisfying | j E | | = .001. SOR was also used to generate this 

l 2 

grid and is in excellent agreement with MGM. The geometry is motivated by 
an impinging jet flow problem that is planned to be run on this grid. The 
flow from the channel interacts with the solid body on the right. The 
computational space has the same shape as the physical space except that 
the body is replaced by a slit. Excluding the channel, the grid consists of 
137 points in the horizontal direction by 97 points in the vertical direction. 
The grid for the channel itself consists of 25 horizontal grid points by 
33 vertical grid points. The slit (and body) are 49 grid points long. 

Corner points on the body and channel have been excluded from the grid. 
Exponential grid spacing was used along various parts of the horizontal and 
vertical boundaries of the grid. In an attempt to preserve this boundary 
spacing in the grid interior non-zero P and Q were used. Although grid lines 
are still bent near the boundaries, they are not bent as much as when P = Q = 0 
was tried. To compute this grid (which had 4 corrective grids, including the 
finest) MGM was "vectorized" on the TI-ASC since it is a parallel processor 
machine. To accomplish vectorization, which cut computing time by a factor of 
six, a four-color relaxation scheme was used (i.e., even points of even rows 
were relaxed simultaneously; odd points of even rows; etc.). With this scheme 
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MGM used 82.781 WU and 45.57 CP seconds to satisfy IJeII = .001 when using 
an overrelaxation factor (RF) of 1.8 on the finest grid and relaxation factors 
of 1.6, 1.4, and 1.2 for the progressively coarser grids. (Varying RF in this 
way improved MGM's performance.) SOR (with the four-color scheme) used 170.0 
WU and 76.68 CP seconds using a relaxation factor of 1.8, which is about 
optimum for this SOR. MGM used 60.641 WU and 36.67 CP seconds to satisfy 
| |E| L = .001 when RF's of 1.6, 1.4, 1.2, 1.0 were used on progressively 
coarser grids (with 1.6 used for the finest grid). With these RF's MGM used 
26.016 WU to satisfy 1 j E j | = .01 compared to 82.0 WU used by SOR with 

RF = 1.8. The parameters 6 = .05, n = .95 were used for MGM which was 
divergent for n < .9. MGM should perform better with a better initial guess 
than used here. (The horizontal straight lines in the initial guess were 
discontinuous at the right-most boundary.) 
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